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I. INTRODUCTION 

Periodic gravitational- wave signals like those originat- 
ing from rotating neutron stars are an important class 
of sources that can be detected by currently operating 
ground-based detectors. Several methods were developed 
to search for such sources and several searches were per- 
formed. This paper continues the series of papers 0- 
Q devoted to studies of data analysis tools and algo- 
rithms needed to perform an all-sky coherent search for 
quasiperiodic gravitational waves. 

The search presented in the current paper is based on 
the maximum-likelihood statistic called the F-statistic 
that we have derived in the Paper I [l| of this series. It is 
known that the coherent search for long observation time 
needed to detect weak gravitational- wave signals from ro- 
tating neutron stars are computationally prohibitive (see 
Q and Paper III of this series [I])- Promising strategies 
are hierarchical semi-coherent methods. In these meth- 
ods data is broken into short segments. In the first stage 
each segment is analyzed using the F-statistic and in the 
second stage the F-statistics from the short segments are 
combined using a certain algorithm. There are several 
methods proposed for the second stage: search for coin- 
cidences among candidates from short duration segments 
@, [1] , stack-slide method H, power flux method |9, E3 > 
Hough transform method 0, [llTtl4| . Recently an opti- 
mal method for the second stage has been found, the 



global correlation coordinate method [la , Il6j | , which ex- 
ploits global parameter space correlations in the coherent 
detection statistic. In our paper we shall present meth- 
ods to optimize the first, coherent stage of a hierarchical 
method. 

The techniques presented in this paper were used in 
the analysis of NAUTILUS bar detector data [ijj and 
are presently used in the analysis of the VIRGO data. 
Alternative techniques for the coherent stage based on 
the F-statistic and their application to the real data can 
be found in Refs. @, H H E3 

The paper is organized as follows. In Sec. [IT] we present 
the noise-free response of a ground-based detector to a 
gravitational-wave signal from a rotating neutron star. 
This response was derived and discussed in detail in Pa- 
pers I [l| and IV [j] of our series. In Sec. IIIII we present 
data analysis tools to perform coherent search of the data 
for a gravitational-wave signal given in Sec. [ill I n Sec. 
1111 Al we present the F-statistic that was derived in Pa- 
per I. We limit ourselves to the case when the observa- 
tion time is an integer multiple of one sidereal day. This 
simplifies some general formulas considerably. In Sec. 
IIII Bl we introduce a simplified approximate model for 
a periodic gravitational-wave signal. This approximate 
signal has the constant amplitude and its phase is pa- 
rameterized in such a way that it is a linear function of 
the parameters. For such a signal the Fisher matrix is 
constant and consequently it is independent of the val- 
ues of the signal's parameters. In Sec. IIII CI we briefly 
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review calculation of the false alarm probability. Section 
ITVl is devoted to construction of the grid of templates in 
the parameter space. The grid solves a certain covering 
problem with a constraint. Our constraint is that the 
nodes of the grid coincide with the Fourier frequencies. 
This allows to use the fast Fourier transform (FFT) al- 
gorithm to compute the ^-statistic at grid nodes, what 
greatly accelerates the calculation. In Sec.[V]we describe 
our package Top2Bary that is used to calculate the posi- 
tion and the velocity of the detector located on the Earth 
with respect to the solar system bary center. In Sec. IV Al 
we introduce various concepts and definitions used in the 
astrometry and in Sec. IV Bl we describe the content of 
our package which is a set of FORTRAN routines. In Sec. 
IVII we present various approximations that we use in the 
calculation of the ^-statistic in order to speed up compu- 
tations. In Sec. VI A we discuss resampling of the time 
series to the barycenter that we need to perform before 
we can apply the FFT. We develop two algorithms: one 
slow and very accurate and the other fast but less ac- 
curate. We compare the two algorithms using the signal 
from Sec. [II] In Sec. VI B we describe interpolation of the 
FFT in the Fourier domain. This interpolation method 
allows to obtain efficiently an FFT that is twice as fine 
as the FFT of original data. In Sec. VI C we describe the 
Nelder-Mead algorithm that we use to find accurately the 
maximum of the ^-statistic. In Sec. I VIII we perform a 
number of Monte Carlo simulations of the computer code 
where we have implemented the methods and algorithms 
from Sees. Ill- VI. In our simulations we investigate how 
well we estimate the parameters of the signal in compar- 
ison to the Cramer-Rao bound. 



II. RESPONSE OF A DETECTOR TO 
A PERIODIC GRAVITATIONAL WAVE 

The dimensionless noise-free response h of a 
gravitational- wave detector to a weak plane gravitational 
wave in the long wavelength approximation [i.e., when 
the size of the detector is much smaller than the reduced 
wavelength X/(2tt) of the wave] can be written as the 
linear combination of the two independent wave polar- 
izations h + and h x , 



h(t) = F+(t)h+{t) + F x (t)h x (t), 



(2.1) 



where F + and F x are the detector's beam-pattern func- 
tions, which are of the form 

F + (t) =sinC(a(£)cos2^ + &(t)sm2^), (2.2a) 

F x (t) = sin ((b(t) cos 2?p - a(t) sm2ip) . (2.2b) 

The beam-patterns F+ and F x are linear combinations 
of sin 2tp and cos 2ip, where tp is the polarization angle of 
the wave. For interferometric detectors the angle £ is the 
angle between the interferometer arms (usually £ = 90°) 
whereas for the case of bars one has to put C = 90° . The 



functions a(t) and b(t) are amplitude modulation func- 
tions, which depend on the location of the detector on 
the Earth and on the position of the gravitational-wave 
source in the sky (described in the celestial coordinate 
system by the right ascension a and the declination 5 of 
the source) . They are periodic functions of time with the 
period of one sidereal day. Analytic form of the functions 
a(t) and b(t) depends on the type of the detector; for the 
case of bar detectors they are explicitly given in Eqs. 
(All) of Ref. Q, whereas for interferometric detectors 
they can be found in Eqs. (12) and (13) of Ref. [l|. 

We are interested in periodic waves, for which the wave 
polarization functions are of the form 

h + (t) = h 0+ cos(4>(t) + 4>o), (2.3a) 
h x (t) = h 0x sm((f>(t) + (f> a ), (2.3b) 

where ho+ and ho x are constant amplitudes of the two 
polarizations and (f>(t) + (f>o is the phase of the wave (with 
(fto being the initial phase of the waveform). The ampli- 
tudes ho+ and hg x depend on the physical mechanism 
generating gravitational radiation. E.g., if a neutron star 
is a triaxial ellipsoid rotating around a principal axis with 
frequency /, then these amplitudes are 



ho+ = 2^0(1 + cos 2 l), 
hox = n o cos l, 



(2.4a) 
(2.4b) 



where t is the angle between the star's angular momen- 
tum vector and the direction from the star to the Earth, 
and the amplitude ho is given by 



16ir 2 G elf 2 



(2.5) 



Here I is the star's moment of inertia with respect to the 
rotation axis, r is the distance to the star, and e is the 
star's ellipticity defined by e = \Ii — h]/!, where I\ and 
I2 are moments of inertia with respect to the principal 
axes orthogonal to the rotation axis. 

We further assume that the gravitational waveform 
given by Eqs. (|2 . 1 1) — (|2 .3[) is almost monochromatic 
around some angular frequency wo, which we define as 
instantaneous angular frequency evaluated at the solar 
system barycenter (SSB) at t = 0. The phase modu- 
lation function <j){t) for such waveform is approximately 
given by 



fc+i 



n • r d (i) 



fc=0 



(fc + 1)! 



5>-, (2.6) 



k=0 



where ujk (k — 1, 2, . . . , s) is the fcth time derivative of the 
instantaneous angular frequency at the SSB evaluated at 
t = 0, no is the constant unit vector in the direction of 
the star in the SSB reference frame (it depends on the 
right ascension a and the declination 6 of the source), 
and rd is the vector joining the SSB with the detector. 
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Approximations that lead to Eq. (|2.6I) are discussed in 
detail in Sec. II B and Appendix A of Paper I. 

Let us associate the following coordinate system with 
the SSB reference frame. The x axis of the system is 
parallel to the x axis of the celestial coordinate system, 1 
the z axis is perpendicular to the ecliptic and has di- 
rection of the orbital angular momentum vector of the 
Earth. In this SSB coordinate system the vector no has 
the components 



n 



1 

cos e sin e 
— sin e cos e 



cos a cos 5 
sin a cos 5 
sin S 



(2.7) 



where e is the obliquity of the ecliptic. The position 
vector r d of the detector with respect to the SSB has the 
components 



rd = 



-""ES 

V o 



1 

cos e sin e 
— sin £ cos £ 




(2- 



where (R%$, Res,0) are the components of the vector 
joining the SSB with the center of the Earth in the SSB 
coordinate system, and {R%, i? E , are the components 
of the vector joining the center of the Earth and the de- 
tector's location in the celestial coordinate system. Mak- 
ing use of Eqs. (|2.7p and (|2.8p one can obtain the explicit 
formula for the scalar product no • r^(t): 

n • r d (t) = cos a cos 5 (i? ES (t) + i?| (t)) 

+ (sin a cos S cos £ + sin 5 sin e) i? Eg (t ) 

+ sin a cos <Si? E (i) + sin<5 i? E (t). (2.9) 



The phase of the waveform given by Eq. (|2.6[) de- 
pends on the angular frequency ujq, the s spin-down 
parameters LUk (k = l,...,s), and on the angles a, 
S (through the vector no). We call the parameters 
(wo, wi, . . . , oj s , a, S) the intrinsic parameters and the re- 
maining ones (/io+, hox , </>0j "0) the extrinsic (or ampli- 
tude) parameters. As we shall see in the following sec- 
tion we only need to search for signals over the intrinsic 
parameter space. The whole signal h depends on s + 7 un- 
known parameters: (/io+, ^ox , </>o 5 "0j ot, <5, u>q,uji, . . . , lu s ). 

The response function h depends on the position of the 
detector with respect to the SSB. This position can be 
determined with a great accuracy using JPL Planetary 
and Lunar Ephemerides DE405/LE405 as described in 
Sec. El The dominant term in the phase 0(t) is u>o t; 
typical gravitational- wave frequency fo :— 2t:/ujo is con- 
tained in the range from a few Hz to a few kHz. The 



1 In the celestial coordinate system the z axis coincides with the 
Earth's rotation axis and points toward the North pole, the x 
and y axes lie in the Earth's equatorial plane with the x axis 
pointing toward the vernal point. 



gravitational- wave signal from a rotating neutron star is a 
nearly periodic signal that is weakly amplitude and phase 
modulated due to the intrinsic variation of star's rotation 
frequency and the motion of the detector with respect to 
the star. Moreover the amplitude of this signal is ex- 
pected to be very small. Consequently detection of the 
signal requires observation time T that is very long with 
respect to the gravitational- wave period Pq := 2tt /ujq. 

Combining Eqs. (|2. 1|) — (|2.3[) together one can decom- 
pose the response h into linear combination of four time- 
dependent components : 



h(t) = ^Ai hi(t), 



(2.10) 



where the functions hi (i = 1, . . . , 4) are of the form 
h^t) = a(t)cos<f>(t), h 2 {t) = b(t) cos (f>(t), 



h 3 (t) = a{t) sin <p(t) , hi (t) = b(t) sin <j>(t) , 



(2.11) 



and the four constant amplitudes Ai (i = 1, ... ,4) are 
given by 

Ai =Aism(, i = l,...,4, (2.12) 

where 

Ax — h 0+ cos 2i[i cos </>o — /iox sin 2-0 sin0 o , (2.13a) 

A~2 — h + sin 2ip cos O + h 0x 008 2^8^00, (2.13b) 

A 3 = — h 0+ cos 2-0 sin O — h ax sin 2ip cos (f>o, (2.13c) 

A4 = — h 0+ sin 2^ sin 0o + ^-ox cos 2ip cos 0o . (2.13d) 



One can invert Eqs. (|2.13[) to obtain formulas for the 
parameters /io+> ^ox? 0Oi and ?/> as functions of the am- 
plitudes Ai. Let us introduce quantities 

A := A\ + Al + A 2 3 + Al, (2.14a) 

D := AiAi -A 2 A 3 . (2.14b) 

Then the amplitudes ho + and hg x can be uniquely deter- 
mined from the relations (we assume here, without loss 
of generality, that ho+ > 0) 



^04 



ho> 



-(a+ \/A 2 -4L> 2 



ign(D)J±(A-VA 2 -4D> 



(2.15a) 



(2.15b) 



The initial phase 0o and the polarization angle ip can be 
obtained from the following equations: 



tan 20o 



2(A 1 A 3 + A 2 A 4 ) 

A3 ' A i A A2 

. ., 2(A x A 2 + A 3 Ai) 
tan 4^ = 

A l + A3 A2 ^4 



(2.16a) 
(2.16b) 



4 



Also Eqs. (|2.4I) can be solved for the amplitude ho and 
the angle i. The result is 



l = arccos(/iox/^o)- 



(2.17a) 
(2.17b) 



In the special case when the star's angular momentum 
vector lies along the line of sight, cost = ±1, and the 
number of independent amplitude parameters is reduced 
to two. In this situation Eqs. (|2.4I) read (upper sign is 
for cos i = +1 and lower sign is for cos t = — 1) 

h 0+ = h 0l h 0x = ±h , (2-18) 

and Eqs. (|2.13p simplify then to 

A 1 = h cos(2ip±<p ), (2.19a) 

A 2 = h sm{2iP±(j)o), (2.19b) 

A 3 = tA 2 , (2.19c) 

A i = ±A 1 . (2.19d) 

III. MAXIMUM-LIKELIHOOD FILTERING 
A. The ^-statistic 



The gravitational- wave signal h given by Eqs. (|2.10[) 
and (|2.1ip will be buried in the noise of a detector. We 
are thus faced with the problem of detecting the signal 
and estimating its parameters. A standard method is 
the method of maximum-likelihood (ML) detection that 
consists of maximizing the likelihood function, which we 
shall denote by A, with respect to the parameters of the 
signal. If the maximum of A exceeds a certain threshold 
calculated from the false alarm probability that we can 
afford, we say that the signal is detected. The values 
of the parameters that maximize A are said to be the 
maximum-likelihood estimators of the parameters of the 
signal. The magnitude of the maximum of A determines 
the probability of detection of the signal. 

We assume that the noise n in the detector is an ad- 
ditive, stationary, Gaussian, and zero-mean continuous 
random process. Then the data x (if the signal h is 
present) can be written as 

x(t)=n(t) + h(t). (3.1) 

The logarithm of the likelihood function has the form 



lnA=(x\h)-±(h\h), 
where the scalar product ( • | • ) is defined by 



(x\y) :=-5R 



x(uj)y*(uj) 



du. 



(3.2) 



(3.3) 



In Eq. (13.3[) tilde denotes the Fourier transform, asterisk 
means complex conjugation, Sh is the one-sided spectral 
density of the detector's noise, and 5i denotes the real 
part of a complex expression. 

We further assume that over the frequency bandwidth 
of the signal h the spectral density Sh is nearly constant 
and equal to So = Sh{^o), where ujq is the frequency of 
the signal measured at the SSB at t = 0. Then the scalar 
products entering Eq. (I3.2|) can be approximated by 



(3.4a) 



(x\h)n— / x(t)h(t)dt, 

JO Jo 

»)4 [ T °(h(t)) 2 dt, 

<-50 JO 



(3.4b) 



where T D is the observation time, and the observation 
interval is (0,T o ). It is useful to introduce the following 
notation 



(3.5) 



x) :=i / ° x{t)dt. 

1 o Jo 



After applying this notation and making use of Eqs. 
the log likelihood ratio from Eq. (|3.2[) can be written as 



(3.6) 



In Sec. Ill of Paper III we have analyzed in detail the 
likelihood ratio for the general case of a signal consisting 
of several narrow-band components. Here we only sum- 
marize the results of Paper III and adapt them to the 
case of our signal (|2.10l) . The signal h depends linearly 
on four amplitudes Ai. The likelihood equations for the 
ML estimators Ai of the amplitudes Ai are given by 



<91nA 
~dA~ 



0, 



1. 



,4. 



(3.7) 



One can easily find the explicit analytic solution to Eqs. 
p.7[) . To simplify formulas we assume that the observa- 
tion time T a is an integer multiple of one sidereal day, 
i.e., T Q — n(2w/Q r ) for some positive integer n, where fi r 
is the rotational angular velocity of the Earth. Then the 
time average of the product of the functions a and b [see 
Eqs. (|2.2I) ] vanishes, (ab) — 0, and the analytic formulas 
for the ML estimators of the amplitudes are given by 



At 



Ao 



(xhi) 



(xh 3 ) 



An 



Aa 



(xh 2 ) 
(b 2 ) ' 
(xhi) 



(3i 



Explicit formulas for the time averages (a 2 ) and (6 2 ) can 
be found in Appendix B of Paper IV. 

The reduced log likelihood function T or the T- 
statistic is the log likelihood function (|3.6[) with the am- 
plitude parameters Ai replaced by their estimators A4. 
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By virtue of Eqs. (|3 . 8[) from Eq. (|3.6[) one gets 



(62) 



where 



Ft 



b • = 



a;(t) a(t) exp[— i(j)(t)} dt, 



x(t) b{t) exp[-i0(i)] di. 



(3.9) 

(3.10a) 
(3.10b) 



The ML estimators of the signal's parameters are ob- 
tained in two steps. Firstly, the estimators of the fre- 
quency, the spin-down parameters, and the angles a and 
S are obtained by maximizing the functional J- with re- 
spect to these parameters. Secondly, the estimators of 
the amplitudes Ai are calculated from the analytic for- 
mulas (|3.8j) with the correlations (xhi) evaluated for the 
values of the parameters obtained in the first step. 



B. A linear model 

In this subsection we introduce a useful approximate 
model of the gravitational-wave signal from a rotating 
neutron star. The model relies on (i) neglecting all spin 
downs in the phase modulation due to motion of the de- 
tector with respect to the SSB; and (ii) discarding this 
component of the vector rd (connecting the SSB and the 
detector) which is perpendicular to the ecliptic. These 
approximations lead to the following phase of the signal: 



Mt) = E 



k=0 



t k+1 



+ aifii(t) + a 2 p 2 (t), (3.11) 



where a% and a 2 are new constant parameters, 

a% := ujq (sin a cos 5 cose + sin 6 sine), (3.12a) 
a 2 := uj cosacosJ, (3.12b) 

and where p,i(t) and p 2 (t) are known functions of time, 
Hl(t) := -(-R|s(*) C0S£ )> (3.13a) 

:=~(iZ|s(t)+i2|(i)). (3.13b) 

We also neglect the slowly varying modulation of the sig- 
nal's amplitude, so finally we approximate the whole sig- 
nal h(t) by 



h(t) = A COS (<fa n (t) + 0o ) , 



(3.14) 



where Aq and 4>o are the constant amplitude and initial 
phase, respectively. The above signal model is called lin- 
ear because it has the property that its phase (|3.11l) is a 
linear function of the parameters. 



It is convenient to represent the linear model of the 
gravitational-wave signal in the following form 



A I 



h(t; 0) = A cos ( E £fc m fc(*) + O 



fc=0 



(3.15) 



where the vector 6 collects all the signal's parameters, 
:= (Aq, 0o i C): with the vector £ comprising the param- 
eters of the signal's phase, £ := (wq, . . . ,u s , a\, 1x2), 
so £ fe = uj k for k = 0, l,...,s, £ s+i = a x , £, s+2 = a 2 ; 
functions mk(t), k = 0, 1, . . . , s + 2, are known functions 
of time t: m k (t) := t k+1 /{k + 1)! for k = 0,1,..., s, 
m s+ i(t) := (ii(t), and m s+2 (t) := p 2 (t)\ finally, M := 
s + 2. 

For the signal f|3 . 1 5[) we will compute the optimal 
signal-to-noise ratio p, 



p:=JW), (3-16) 
and the components of the Fisher information matrix T, 

— ( — 



ki 



\B9 k 



dh \ 



(3.17) 



It is reasonable to assume that the observation time T Q 
is much longer than the period Pq = 2tt/ujq of the grav- 
itational wave (typically Pq < 0.1s and T Q > lday). As 
a consequence 



(cos[n0un(t)]) « 0, (sin[n0i in (t)]) « 0, 



(3.18) 



for any positive integer n. Making use of these approxi- 
mations one easily computes from Eq. (|3.16[) the signal- 
to-noise ratio, 



(3.19) 



and from Eq. (|3.17l) the components of the Fisher infor- 
mation matrix, 



T A A 
^A <p 



, 2 , 1 0000 

^0 



p 2 (mfc) 



0, fc = 0, 
fc = 0,.. 



p 2 (m k m e ) , k,£ = Q,...,M. 



(3.20a) 
(3.20b) 
(3.20c) 
(3.20d) 



Assuming that the signal (|3.15p is buried in the sta- 
tionary and Gaussian noise, one easily computes its T- 
statistic, 



SqT q 

where x(t) are the data 



M 



x(t) exp I - i E £fc m ft (t) ) dt 

(3.21 
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C. False alarm probability 

Let us calculate the autocovariance function C of the 
F-statistic (|3.21[) in the case when data is only noise. It 
is defined as 

C(£,0 := E [^(C)-F(C')] - E [.F(0]Eo[.F(OL (3-22) 

where Eo is the expectation value when data is only noise. 
We find that E [F(£)] = 1 and that 



C(t) « ^cos y^T k m k (tyj \ + / sin Tfcmfe \ , 

(3.23) 

where r := £ — £ . Thus the autocovariance function 
depends only on the difference of the parameters at two 
points but not on the parameters themselves. 

The autocovariance C attains its maximum value equal 
1 when t = 0. Let us consider Taylor expansion of C 
around the maximum up to terms quadratic in r, 



C(t)«1 + £ 



8C(t) 



T=0 



It. 



3 2 C{t) 



2 ^ dT k dri 



T = 



As C attains its maximum for r = 0, we have 
dC(r) 



0. 



(3.24) 



(3.25) 



Let us introduce the symmetric matrix G with elements 



Ghi :— — — 



1 8 2 C(t) 



2 drkdn 



(3.26) 



T = 



One can show that G = Y. where Y is the reduced Fisher 
matrix defined by 



fki := (m k mi) - (m k ) (mi) 



(3.27) 



For the linear phase model the components of the reduced 
Fisher matrix are constants independent of the values of 
the parameters. Making use of Eqs. (|3.25l) - (|3.27j) . the 
Taylor expansion (|3.24p can be written in the form 



C(t) « 1 - y^ffcl TfcTl. 



(3.28) 



We define now the correlation hypersurface of the 
statistic F by the requirement that the autocovariance 
C attains some constant value Cq on it: 

C(t)=C . (3.29) 

This equality, by virtue of Eq. (|3.28[) , can be written as 

I- C . (3.30) 



T k Tit 

k,£ 



Equation (|3.30[) defines an M-dimensional hyperellipsoid. 



The main idea is to divide the space of the phase pa- 
rameters £ into elementary cells which boundary is de- 
termined by Eq. p.30[) . We choose the value Co = 1/2. 
We estimate the number N c of elementary cells by divid- 
ing the total Euclidean volume Vtotai of the parameter 
space by the Euclidean volume V^ c n of the correlation 
hyperellipsoid, i.e., we have 



Nr. 



Vt 
V, 



total 



cell 



(3.31) 



where the Euclidean volume of one elementary cell equals 



Vceli 



T M/2 



T(M/2 + l)VdetG' 



(3.32) 



here Y denotes the Gamma function. 

The values of the statistic T in different cells can be 
considered as independent random variables. We approx- 
imate the probability distribution of T in each cell by the 
probability distribution po (J 7 ) of T when the signal is ab- 
sent. When the signal is absent, 2 J 7 has a \ 2 distribution 
with four degrees of freedom. The false alarm probabil- 
ity Pp for a given cell is the probability that T exceeds 
a certain threshold T$ when there is no signal; for \ 2 
distribution with four degrees of freedom we have 



P F (T ) = (1 + Jo) exp(- Jo) 



(3.33) 



The probability that T does not exceed the threshold 
Tq in a given cell is 1 — Pp(Jo), where Pf(Fq) is given 
by Eq. (|3.33p . Consequently the probability that T does 
not exceed the threshold J-q in all the N c cells is (l — 

Pp(J r o)) Nc . The probability a that J- exceeds T§ in one 
or more cell is thus given by 



1- (l-Pp(Jb)) 



(3.34) 



This is the desired false alarm probability. Inverting the 
formula (|3.34[) we can calculate the threshold value 
corresponding to a chosen false alarm probability a. The 
expected number of false alarms Np is given by 



Np = N c P F (T ). 



(3.35) 



IV. GRID IN THE PARAMETER SPACE 

In order to search for a signal in the noise of the detec- 
tor we need to construct a grid in the space of the signal's 
parameters. We define a grid in such a way that for any 
possible signal there exists a grid point in the parameter 
space such that the expectation value of the .F-statistic 
for the parameters of this grid point is greater than a 
certain value. 

In the construction of the grid we employ the ap- 
proximate linear model of the signal introduced in Sec. 
IIII Bl The F-statistic for this signal is given in Eq. 
(|3.21j) . The expectation value of the F-statistic, when 
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the signal is present in the data [i.e., when the data 
x(t) = n(t) + h(t;9), where 6 = (A ,(j> ,£) collects the 
signal's parameters], is equal 

Ea [F[n{t) + h(t; 6);$']] = 1 + ^C(r), (4.1) 

where p is the optimal signal-to-noise ratio computed in 
Eq. (|3.19j) and C is the autocovariance function given 
in Eq. (|3.23[) ; the vector r := £ — where are the 
phase parameters of the template. The function C has 
the maximum equal to 1 for t — 0. The autocovariance 
function C is equal to the square of the match function 
M originally defined by Owen (20| . 

To construct the grid we first choose the minimum 
value of the correlation that we can accept. We denote 
this value by C Q . (Let us note that C = MM\ where MM 
is the minimal match introduced by Owen [20j.) Then we 
introduce, as in Sec. IIII C[ the correlation hypersurface 
of the statistic T by the equality C(r) = Co, which, after 
making the Taylor expansion of C up to the second order 
terms in r, is described by Eq. (|3.30p . 



A. The covering problem with constraints 

The problem of constructing a grid in the parame- 
ter space is equivalent to the so called covering problem 
[22l [23| : we want to cover (M + l)-dimensional param- 
eter space with identical hyperellipsoids (|3.30[) in such a 
way that any point of the space belongs to at least one 
ellipsoid. Moreover, we look for an optimal covering, i.e., 
the one having smallest possible number of grid points 
per unit volume. The covering thickness 9 is defined as 
the average number of ellipsoids that contain a point in 
the space. The optimal covering would have minimal 
possible thickness. 

Let us introduce in the parameter space the new set of 
coordinates x = (xq, . . . , xm), defined by the equality 

t = Mx, (4.2) 

where the transformation matrix M is given by 

M = U D - 1 . (4.3) 

Here Do is the diagonal matrix whose diagonal compo- 
nents are square roots of eigenvalues of L, and Uo is a 
matrix whose columns are eigenvectors of L, normalized 
to unity. One can show that Uo is an orthogonal matrix, 
U^ 1 = Up (superscript 'T' denotes matrix transposition) 
and 

f = U DgUj. (4.4) 

Hyperellipsoid (|3.30l) in coordinates x reduces to the 
(M + l)-dimensional sphere of radius R := y/1 — Co- 
Therefore, the optimal grid can be expressed by means 
of a sphere covering. 



In general, the thinnest possible coverings are known 
only in dimensions 1 and 2. In dimensions up to 5 the 
thinnest lattice coverings are known [see Eq. f|4.7[) below 
for the definition of a lattice], while in many higher di- 
mensions the thinnest known coverings are lattices (22j. 
From this point on we consider only lattice coverings. 
However, these general results cannot be easily adopted 
to our case. For computational reasons, we would like 
the nodes of the grid to coincide with Fourier frequen- 
cies, so that we can use the FFT algorithm to calculate 
the .F-satistic efficiently. 

Our grid should meet the following constraint: one of 
its basis vectors needs to lie on the frequency axis and 
have given length. In other words, we look for the optimal 
covering with one of the lattice vectors fixed. We denote 
this vector by 

a o = (A Po ,0,...,0), (4.5) 

where Apo is the fixed frequency resolution of our proce- 
dure. There is another constraint to be met in an all-sky 
search: 

a t = (aio,aa,...,a is ,0, 0), i = l,...,s. (4.6) 

Having this constraint satisfied, we greatly reduce the 
computational overhead of resampling the data to the 
barycentric time (see Sec. I VI Aj) . 

As far as we know, the general solution to the covering 
problem with constraints is not known. Starting from the 
hypercubic covering (i.e., having all the lattice vectors or- 
thogonal), a covering satisfying both constraints can be 
constructed, for the signal (|3.15p with any number of pa- 
rameters (see pljV However, in higher dimensions it may 
be several times thicker than the thinnest unconstrained 
lattice known. An improved construction is proposed 
here, which takes as a starting point the thinnest lattice 
covering known in a given dimension, and applies a se- 
quence of modifications to satisfy the constraints. We 
will refer to this lattice as the optimal covering. 

B. Optimal lattice 

A lattice can be conveniently defined as a set of all 
linear combinations of its basis vectors with integer 
coefficients: 

A= j^fc.a, : h G zj . (4.7) 

Given lattice A, its fundamental parallelotope is the set 
of all points of the form J2i 9i a i, with < 0i < 1. Fun- 
damental parallelotope is one example of elementary cell. 
The thickness 9 of a lattice covering is equal to the ra- 
tio of the volume of one hyperellipsoid to the volume of 
fundamental parallelotope. 

For any lattice point G A, the Voronoi cell around 
P, is defined as 

V(Pi) = {t : C(t - Pi) > C(t - P,-) for all 1} , 

(4.8) 
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where C(r) is the Taylor expansion f|3 . 28[) of the autoco- 
variance function. All Voronoi cells of any lattice A are 
congruent, disjoint, and their union is the whole space. 
Voronoi cell is another example of elementary cell and 
is sometimes called Wigner-Seitz cell or Brillouin zone. 
The Voronoi cell of A is inscribed into the correlation 
ellipsoid (j33Ul) . 

Let A be any lattice with basis vectors (ao , ai , . . . ). 
The square of minimal match of A is 

MM 2 (A) = inf CYt-PA (4.9) 

where P, can be any lattice point. Let £ g V(Pi) be the 
value for which the minimum in (j4.9|l is found. The func- 
tion C{t — Pi) has at the point £ its absolute minimum 
inside the Voronoi cell V(Pj), and C is a deep hole of A. 
Note that the deep hole must be one of the vertices of 



the Voronoi cell. It makes Voronoi cells especially useful 
for calculating minimal match of a given lattice. 

We can now outline the construction of an optimal cov- 
ering in the parameter space. Given the value of Co, we 
look for the thinnest possible lattice covering A, satisfy- 
ing 

MM 2 (A) = C* . (4.10) 

As a starting point, we consider the thinnest lattice cov- 
ering known. It is determined by the number of phase 
parameters. For example, the thinnest covering of 4- 
dimensional space is the so called Voronoi's principal 
lattice of the first type A\ [22| . having the thickness 
$min = 1-7655. The generator matrix (a matrix whose 
rows are the basis vectors) of this lattice reads 
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(4.11) 



2V5/ 



where R is the covering radius. 

Let 1 be any lattice vector of such that |1| > 
|M~ 1 a |, where M and a are given by Eqs. (|4.3p and 
(|4.5p . respectively. The vector 1 is a linear combination 
of rows of (|4.1ip . In order to construct an optimal lattice, 
satisfying constraint (|4.5|) . we perform two scaling oper- 
ations on A\ : shrinking in the direction of 1 by the factor 
/i = | M ~ 1 ao I / 1 1 1 < 1 and expanding in all directions per- 
pendicular to 1 by the factor v > 1. The generator matrix 
of the shrunk lattice is 

Mi = M Oi diag( M , 1, 1, 1) Or 1 , (4.12) 

where Oi is an orthogonal transformation defined by the 
condition that Oi (1,0, 0,0) = 1/|1|. 

The expansion factor v is defined by the condition that 
the covering radius R remains unchanged on scaling. In 
general, its value can be determined numerically using 
the following iteration: 

M i+ i = M l Oidiag(l,^,^,i/. i )Or 1 I for « = 1,2,..., 

(4.13) 

where Vi = R/{Ri sin^), Ri is the covering radius of the 
lattice given by M^, and <pi < ir/2 is the largest angle 
that a half-line starting at the origin and containing a 
deep hole of can form with the direction of 1. The 
above procedure converges after several iterations, and 
the expansion factor is formally v = Yi^Li v %- 

After scaling, the lattice thickens by a factor 1 / (/xi/ 3 ) > 
1, depending on the choice of 1. Note that /iz/ 3 —> as the 



length of 1 increases. By enumerating the lattice vectors 
of A\ in order of increasing magnitude, one can find the 
optimal 1, such that l/(^ 3 ) is minimal. 

The generator matrix of an optimal lattice, satisfying 
the constraint (|4.5[) . is 

M opt = M OiSO2 1 M T , (4.14) 

where O2 is an orthogonal transformation satisfying 
02(1,0,0,0) = e /|e | with e = M _1 a , S is a diago- 
nal matrix with elements The lattice vector 
1 is chosen in such a way that det S is maximal. The 
thickness of this lattice is 



For the case of the Virgo antenna, observational time 
To = 2 days, frequency 750 Hz, effective bandwidth 1 Hz, 
and Co = 3/4, the resulting lattice is thicker than A\ only 
by 20%, 9 ~ 2.1, reducing the number of templates to 
roughly 10 millions (the corresponding hypercubic lattice 
would have 25 millions of templates). 

The constraint (|4.6p can be satisfied without increasing 
the lattice thickness by a 4-dimensional rotation O3 such 
that Oseo = eo and the vector M03020j~ 1 f is orthogonal 
to the {0.1,0.2) plane, where f is a lattice vector of A\, 
f \ 1. The generator matrix of a lattice satisfying both 
constraints is now 

M opt = M Oi S O^T 1 Or 1 M T . (4.16) 
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C. Two-dimensional example 

Let us explain the construction of an optimal lattice 
on a simple, two-dimensional example. We consider here 
the signal 



t 



h(t;A ,<j} ,u} ,u}i) = A cos (^o^r + w i 

(4-17) 

where T a is the observation time and the observation in- 
terval is (0,T o ). The phase of the signal (|4.17p depends 
on two dimensionless parameters uiq, lj\ (and on the ini- 
tial phase parameter (f>o). The reduced Fisher matrix for 
this signal reads 



(4.18) 




The correlation hypersurface (|3.30p is now an ellipse: 

(n 



T2 



12 



— = 1 - C . 
180 



(4.19) 



Correlation ellipse is shown on both parts of Fig. Q] Fre- 
quency resolution when calculating ^-statistic is now 7r, 
therefore we require the first basis vector of a lattice to 
be equal to 



a = (tt,0), 



(4.20) 



in order to satisfy the constraint (14.51) . 

The thinnest possible covering in two dimensions is 
the hexagonal lattice, A^. The Voronoi cell of A\ is a 
regular hexagon, and the covering thickness is f?hcx = 
2tt/(3V3) = 1.2092. The generator matrix of A* 2 is 



1 = W3(i % 

-2 2 



(4.21) 



where R = \Jl — Co is the covering radius. In order to 
cover the parameter space with ellipses (|4.19p . we have 
to transform the hexagonal lattice to the r = {t\,T2) 
coordinates, as given by (|4.2j) . The generator matrix of 
an optimal, hexagonal covering in r coordinates is 



Mhcx — 



(4.22) 



The resulting lattice is shown in the upper part of Fig. 
[TJ along with the correlation ellipse. The Voronoi cell is 
a hexagon inscribed into the correlation ellipse (any lat- 
tice A constitutes a covering if and only if the correlation 
hypersurface completely includes its Voronoi cell). Note 
that lattice points do not coincide with Fourier frequen- 
cies, represented by open circles. 

The requirement (|4.20j) can be satisfied by the means 
of transformation (I4.14j) . at the cost of increasing lat- 
tice thickness. We find that |M^ 1 a | = tt/(2 v / 3). For 
C = 3/4 (then R = 1/2), the best choice of 1 is 
1 = -\/3/4(3, -\/3). The hexagonal lattice in coordi- 
nates Xi is shrunk in the direction of 1 by the factor 
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FIG. 1: Up: Hexagonal covering generated by (|4.22[1 . with 
Co = 3/4 and R — 1/2 (dark points). Fourier frequencies 
are represented by open circles. Down: Constrained lattice, 
satisfying the condition (|4.20p . Any point in the parameter 
space belongs on average to 1.8171 ellipses. 



1 ao|/|l| = 7r/(3v / 3), then expanded by the factor 



M = 

v = Vl08 — 7r 2 /9, and finally rotated in such a way that 
1 coincides with M _1 ao. 

The lattice such obtained is shown in the lower part 
of Fig. Q] It has the thickness 9 = 6h cx /(^v 3 ) = 
18/^108 - 7T 2 S 1.8171. Note that the Voronoi cell has 
changed. It covers more than half of the area of the cor- 
relation ellipse and is inscribed in it. 



V. POSITION AND VELOCITY OF 
A DETECTOR WITH RESPECT TO THE SOLAR 
SYSTEM BARYCENTER 

In this section we describe algorithms and procedures 
needed to compute, for a given time and a geographi- 
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cal location of a detector, the corresponding vectors of 
position and velocity of the detector with respect to the 
Earth barycenter and of this barycenter with respect to 
the solar system barycenter (SSB). The sum of these vec- 
tors represents the position and velocity of the detector 
referred to the SSB. 

We have incorporated the algorithms and procedures 
described in the present section into a set of FORTRAN 
subroutines which we have called the Top2Bary package. 



A. Overview of transformations involved 

1. Reference systems and frames 2 

Reference data for positional astronomy, such as the 
data in barycentric planetary ephemerides, are now spec- 
ified within the International Celestial Reference Sys- 
tem (ICRS). The ICRS is a coordinate system whose 
origin is at the SSB and whose axis directions are ef- 
fectively defined by the adopted coordinates of 212 ex- 
tragalactic radio sources which are assumed to have no 
observable intrinsic angular motions. Thus, the ICRS 
is a "space-fixed" system (more precisely, a kinemati- 
cally non-rotating system) without an associated epoch. 
However, the ICRS closely matches the conventional dy- 
namical system defined by the Earth's mean equator and 
equinox of J2000.0; the alignment difference is at the 0.02 
arcsecond level, negligible for many applications. The list 
of radio source positions that define ICRS for practical 
purposes is called the International Celestial Reference 
Frame (ICRF). 

The position and velocity 3-vectors taken from the JPL 
DE405/LE405 ephemeris are in equatorial rectangular 
coordinates referred to the SSB. The reference frame for 
the DE405 is the ICRF; the alignment onto this frame, 
and therefore onto the ICRS, has an estimated accuracy 
of a few milliarcseconds, at least for the inner-planet 
data. 

The DE405 was developed using T ep h, a barycentric 
coordinate time |3lj . T ep h is rigorously equivalent to 
Barycentric Coordinate Time (TCB) in a mathematical 
sense, differing only in rate: the rate of T ep h matches the 
average rate of TT (Terrestrial Time, or TDT), while the 
rate of TCB is defined by the SI system. The IAU time 
scale, Barycentric Dynamical Time (TDB), often (but er- 
roneously) considered to be the same as T op h, is a quan- 
tity that cannot be physically realized, due to its flawed 
definition. So, in fact, the use of the name TDB actu- 
ally refers to quantities based on or created with T ep h 
(because of this, the IAU Working Group on Nomen- 
clature for Fundamental Astronomy has recommended 
changing the definition of TDB to be consistent with 
that of T ep h). Astronomical constants obtained from 



2 We employ here the recent recommendations of the IAU |2g( . 



ephemerides based on T cp h (or TDB) are not in the SI 
system of units and must therefore be scaled for use with 
TCB or other Si-based time scales. 

The epoch J2000.0 is the epoch of 2000 January 1, 
12 h TT (JD 2451545.0 TT) at the geocenter ("J2000.0 
system" is shorthand for the celestial reference system 
defined by the mean dynamical equator and equinox of 
J2000.0). The coordinate system defined by the "equa- 
tor and equinox of J2000.0", can be thought of as either 
barycentric or geocentric. 

It is also worth noting that the recent IAU resolu- 
tions do not describe the proper reference system of the 
observer — the local, or topocentric, system in which most 
measurements are actually taken. The resolutions as 
adopted apply specifically to Einstein's theory of grav- 
ity, i.e., the general theory of relativity. 



2. Time scales 

It is assumed that the time, associated with obser- 
vational data we are dealing with, is the Coordinated 
Universal Time (UTC) as disseminated by international 
time services. The UTC scale since 1972 is essentially 
uniform, except for occasional 1 second steps (leap sec- 
onds) introduced internationally to compensate for the 
variable Earth rotation. By 2010 there were 24 leap 
seconds. Earlier, 1961 to 1971, the adjustments were 
continuous. UTC devoid of these adjustments is called 
the International Atomic Time, TAI. Addition of the dif- 
ference TAI-UTC to UTC converts it to TAI which is 
uniform. 3 TAI in turn differs from the Terrestrial (Dy- 
namical) Time, TDT or TT (normally used to describe 
celestial phenomena by astronomers), only by a constant 
term: 

TDT = TAI + 32.184 s. 

As already explained, the time argument of the JPL 
positions and velocities of celestial objects is in principle 
T ep h or the barycentric coordinate time. This time scale 
in practice can be equated with the Barycentric Dynam- 
ical Time, TDB (in spite of the noted inadequacy in its 
definition). The TDB differs from the TDT only by small 
periodic terms, which to sufficient accuracy are usually 
simplified to only two largest terms: 

TDB - TDT = (0.001658 sin g + 0.000014 sin 2g) s, 

where g = (357.53 + 0.9856003 (JD-2451545.0)) degrees 
and JD is the Julian Date equal to MJD + 2400000.5, 
MJD being the Modified Julian Date. So essentially, the 



3 The TAI— UTC differences are available in a tabular form at 
hpiers.obspm.fr/eoppc/bul/bulc/UTC-TAI.history The func- 
tion tai_ut of our Top2Bary package reads similar file and returns 
the difference calculated for given UTC Julian Date. 
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two scales do not have to be distinguished and for many 
practical purposes can be assumed equal, TDB=TDT. 4 
One can relate given UTC to barycentric positions of 
all the major celestial bodies of the solar system as ob- 
tainable from the JPL ephemeris. 5 However, to relate 
the position of a point on the Earth to the geocentric 
ICRF (i.e., the same frame as the barycentric ICRF ex- 
cept for the origin of axes which is now at the barycenter 
of the Earth) one has to use yet another time scale — the 
rotational time scale UT1, which is nonuniform and is de- 
termined from astronomical observations. The difference 
UT1— UTC, the value of which is presently maintained 
by international services within ±0.90 s, is taken from 
the International Earth Rotation Service (IERS) tabu- 
lations of daily values publicly available as eopc04 . yy 
files. 6 The UT1 time can be readily converted to the 
Earth rotational angle or the sidereal time. 7 



3. Location coordinates and velocities with respect to the 
Earth barycenter 

To be able to express coordinates of a point on the 
Earth surface in the Earth centered ICRF, it is necessary 
to know orientation of the Earth in space. The primary 
effects that should be taken into account are: diurnal 
(variable) rotation, precession and nutation of the Earth 
rotational axis, and polar motion. 

The precession is accounted for by applying standard 
astronomical theories. We use the new IAU theory [29| . 
The nutation also could be computed basing on a theory, 
but the DE405 has it in the numerical form, so we just 
read the nutation angles, A?/> and Ae. These nutation 
angles are not the same as defined in the newest IAU nu- 
tation theory, so when highest precision is required the 
celestial pole offsets, dip and de, must also be added to 
these angles (in the past the magnitude of these offsets re- 
mained below 0.1 arcsecond). For past years, since 1962, 
these two offsets are included in the already mentioned 



4 In the Top2Bary package we do distinguish them, but the men- 
tioned possibility can be made effective by setting the iTDB op- 
tion, in the useTop2B.cfg configuration file, to 0. 

5 These relations are incorporated into the Top2Bpv subroutine, 
which also calls the tai_ut function that reads a file which con- 
tains details of the UTC adjustments. 

6 Here yy stands for a two-digit year number (e.g., 99 for 1999, and 
06 for 2006). The eopc04.yy files are available at the IERS Inter- 
net site www. iers . org/MainDisp . csl?pid=36-25788&prodid=22 
In the Top2Bary package the UTC to UT1 conversion takes place 
in the sitePV subroutine, which calls a polar motion routine, 
polmot. The latter returns the UT1— UTC difference interpo- 
lated (between two nearest midnight values read from appropri- 
ate eopc04 file) to the UTC given, along with other parameters 
of Earth axis motion (the terrestrial and celestial pole offsets) 
similarly interpolated. 

7 In the Top2Bary package this conversion is performed in the sid 
function, which is a function of UT1 and the location geograph- 
ical longitude corrected for the polar motion. 



eopc04 . yy files. 8 

The remaining two effects, the Earth variable rotation 
and polar motion, are unpredictable for a remoter future, 
so also observational data must be used. The data nec- 
essary for reduction are taken from the eopc04 . yy files 
as well. The polar motion can be taken into account by 
modifying the conventional geographical coordinates of 
a point on the Earth [see Eqs. (5.1) in |4J] or modify- 
ing the rectangular coordinates corresponding to these 
conventional geographical coordinates: 

x = r a cos A , (5.1a) 
y = r Q sin A D , (5.1b) 
z = bsmty + /i sin 0o, (5.1c) 

where <fi is the conventional geographical latitude, A 
— the conventional longitude, h — the height above the 
Earth ellipsoid, *f? = arctan(&tan0 o /et) — the reduced 
latitude, r = a cos ^ + h cos <f> is the equatorial com- 
ponent of the radius vector, and a = 6378.140 km and 
b = a(l — 1/f) are the semiaxes of the ellipsoid (the 
flattening / is taken equal to 1/0.00335281, which is the 
NOVAS value 9 ). 

In this version we have adopted the second possibility 
(i.e. correcting of the Cartesian coordinates) using the 
following relations: 

x — Xq P^Zqj (5.2a) 
y = y +P y z 0l (5.2b) 
z = z + P x x a - P y y , (5.2c) 

where P x and P y are the IERS coordinates of the pole, 
with respect to the Conventional International Origin (to 
which the 'conventional' geographical coordinates refer), 
converted to radians. 

These (x, y, z) coordinates are expressed in the ter- 
restrial reference frame with the x axis directed toward 
the Greenwich meridian. To relate them to the celestial 
frame, the ICRF, the rotational angle of the Earth must 
be taken into account. This is done through conversion 
of UTC to UT1 (as described in the preceding subsec- 
tion). The UT1 time serves to calculate the apparent 
(or true) local sidereal time (returned by the sid func- 
tion), which includes the nutational component (nutation 
in longitude corrected for the corresponding IERS celes- 
tial pole offset also taken from the IERS eopc04 files). 



8 Our package normally adds these offsets, but the user may 
change this option by setting the NutSid parameter in the config- 
uration file useTop2B . cf g. Note, however, that doing so he will 
affect not only nutational transformations but also computing of 
the equation of equinoxes which enters formulas for the sidereal 
time and thus affects the rotation angle of the Earth. 

9 The user may change the value of the flattening / to the IAU 
value of 298.257, by changing the NovF parameter in the config- 
uration file. 
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The apparent local time is advanced with respect to the 
mean Greenwich sidereal time by the true location lon- 
gitude [calculated as A = arctan(y/cc) with proper choice 
of one of the four quadrants] and the mentioned muta- 
tional component. This component is equal to the so 
called equation of equinoxes (Aip + dip) cose plus a very 
small correction to this equation (which amounts to less 
than 0.003" and depends on the mean longitude of the 
ascending node of the Moon). 10 

So computed local sidereal time 9 is finally used to find 
rectangular coordinates of the location in the geocentric 
ICRF: 

X = r e cos 6, (5.3a) 
Y = r e sin0, (5.3b) 
Z = z, (5.3c) 

where r e = \Jx 2 + y 2 is the equatorial component. 

At this point the location velocity due to Earth rota- 
tion can be computed. The location motion relative to 
the Earth barycenter is represented by a vector of con- 
stant length [in principle, v — 2irr e / (sidereal day) = 
fi r r e , where f2 r is the Earth angular rotation speed] and 
directed always towards the east in the topocentric refer- 
ence frame. This diurnal velocity vector has the following 



Cartesian components: 

V x = v cos(9 + tt/2) = — v sin (9, (5.4a) 

V y = v sin(# + tt/2) = +v cos 9, (5.4b) 

V z = 0, (5.4c) 



where the numerical value of v Q is calculated with the 
NOVAS constant of the Earth angular rotation speed, 
fl T = 2?r/(sidereal day in TAI seconds)= 7.2921151467 x 
10" 5 rad/s. 11 

Since our approach is essentially classical, these Carte- 
sian coordinates (X,Y,Z) and velocities (V x ,V y ,V z ) are 
naturally referred to the frame of equator and equinox of 
date. Therefore they are further nutated and precessed 
(in this order) back to the standard epoch J2000.0. 12 



tal solar system ephemerides from the Jet Propulsion 
Laboratory (JPL). The latest JPL Planetary and Lunar 
Ephemerides, DE405/LE405 or just DE405, were created 
in 1997 and are described in detail in Ref. [3l|. 13 The 
DE405 ephemeris is based upon the ICRF. 14 It consti- 
tutes of a set of Chebyshev polynomials fit with full preci- 
sion to a numerical integration over 1600 AD to 2200 AD. 
The JPL package allows to obtain the rectangular coordi- 
nates of the Sun, Moon, and nine major planets anywhere 
between JED (i.e., Julian Ephemeris Date) 2305424.50 
(1599 Dec 09) and JED 2525008.50 (2201 Feb 20). Be- 
sides coordinates, it includes nutations and librations. 15 
The ephemeris gives separately the position and veloc- 
ity of the Earth-Moon barycenter and the Moon's posi- 
tion and velocity relative to this barycenter. The Earth 
position and velocity vectors (relative to the Earth-Moon 
barycenter) are thus calculated as a fraction (involving 
the masses of the two bodies) of the Moon's vectors and 
opposite to them. For example, the x coordinate of the 
Earth barycenter with respect to the SSB is obtained as: 

xe — sem — xm/82. 30056, 

where 'EM' and 'M' subscripts refer to the Earth-Moon 
barycenter and the Moon, respectively, and the numerical 
denominator is equal to the Earth-plus-Moon to Moon 
ratio of masses taken from the JPL ephemeris. Similar 
expressions pertain to the t/e and ze coordinates. Since 
the DE405/LE405 coordinates are given in the J2000.0 
reference frame, the position and velocity of the Earth 
barycenter so obtained need not be nutated nor pre- 
cessed. 16 

Finally, the velocity vector of the motion of the Sun 
towards its apex (with the speed of 20 km/s) can be 
optionally added to the Earth barycenter velocity. The 
direction of solar apex is assumed at 18 h in right ascen- 
sion and 30° in declination in the frame of equator and 
equinox of J1900.0. Therefore this direction has to be 
precessed from that epoch to J2000.0. 17 



4- Bary centric coordinates of the Earth (JPL ephemeris) 

For computing the coordinates of the Earth barycen- 
ter, relative to the SSB, use is made of the fundamen- 



Our configuration file allows the user to neglect this small cor- 
rection or even use the mean sidereal time instead. 

This constant corresponds to NOmega (a parameter listed in the 
configuration file) set to 1, and can be changed to the IERS 
value of 7.292115 X 10" 5 (NOmega = 0) or to 2tt/(24 x 3600) X 
1.002737909350795 = 7.2921158553 X 10~ 5 (NDmega = 2). All 
the above conversions are done in the sitePV subroutine. 
This is performed in the Top2Bpv subroutine by calling the RemNut 
and prexyz procedures, separately for the position vector and 
velocity vector. 



They are available via the Internet 
( ssd. jpl .nasa. gov/?planet_eph_export ) or on 

CDrom (from the publisher: Willmann-Bell, Inc.; 
www. willbell . com/software/ jpl .htm ). 

An earlier popular ephemeris DE200, which has been the basis of 
the Astronomical Almanac since 1984, is within 0.01 arcseconds 
of the ICRF. 

Our routines do make use of the JPL nutation in longitude and 
in obliquity after correction for (addition of) the IERS celestial 
pole offsets. We have used only a 21-year (1990 to 2010) subset 
of the original ephemeris. 

The position and velocity vectors of the Earth are computed by 
calling the EarthPV subroutine, which in turn calls the original 
JPL STATE routine slightly modified for our needs. This subrou- 
tine reads the JPL planetary ephemeris file named DE405'90. '10 
and interpolates the data to the specified epoch. 

Since sky positions in astronomical catalogues are not corrected 
for the apex motion, this component is actually only optionally 
included in the Top2Bary package (corresponding iSunV option is 
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polmot 



VI. SEARCH ALGORITHM 

In the case of all-sky searches for gravitational-wave 
signals from rotating neutron stars the parameter space is 
very large and it is important to calculate the ^-statistic 
as efficiently as possible. 



INTERP 



eps 



DATA 



A. Resampling 



FIG. 2: Block diagram of the Top2Bary module. 



B. Structure of the Top2Bary module 



The above described algorithms and procedures were 
implemented in a module named Top2Bary, consist- 
ing of about 900 lines of FORTRAN code. The 
overall structure of the module is shown in Fig. [2] 
The Top2Bary module contains the following sub- 
programs: Top2Bpv — main subroutine (to be called 
by user); init — sets some constants and parame- 
ters; sitePV — computes location geocentric vectors; 
EarthPV — computes Earth barycentric vectors; STATE — 
modified JPL STATE that reads DE405; INTERP— original 
JPL subprogram; prexyz — precesses rectangular coordi- 
nates using PREnew; PREnew — computes general preces- 
sion (within the new IAU theory 18 ); RemNut — eliminates 
nutation from a vector; nutatJPL — returns JPL nutation 
angles and mean obliquity; eps — calculates mean obliq- 
uity; sid — computes local sidereal time (mean or ap- 
parent); tai_ut — calculates the difference of TAI— UTC; 
polmot — interpolates IERS UT1— UTC and pole offsets; 
DATA — finds Gregorian or Julian calendar date from JD. 

The module expects the presence of the following 
data files: DE405 ' 90 . ' 10 — JPL DE405 ephemeris binary 
file spanning the years 1990 to 2010; tai-utc.dat — 
table of the TAI-UTC differences (ASCII); eopc04.yy— 
IERS files, one per yy-year, for desired years (ASCII); 
useTop2B. cfg — optional configuration file (ASCII). 19 



set to in the code). If desired, the iSunV option can be changed 
by the user just by setting a nonzero value in the configuration 
file, useTop2B . cfg. In this case the velocity vector towards the 
apex is added to the Earth velocity vector (position remaining 
unaffected). 

18 The once adequate word 'new' is misleading in view of recent rev- 
olutionary changes of concepts since by this name really referred 
is here that old IAU 1976 theory 

19 Except for the last file which is looked for in the current direc- 
tory, all the other files must be placed in the /Top2Bary/EphData 
subdirectory. All the ASCII files must be formatted in Windows 
style (i.e., the lines must be terminated with CR/LF) and not 
UNIX style (with LFs only). This is important for a user who up- 
dates or modifies existing files or downloads new eopc04 . yy files 
from the IERS site, where they are stored in the UNIX format. 



The detection statistic T of Eq. Q3.9p involves integrals 
given by Eqs. (|3.10p . Let us consider the integral (I3.10a|) 
[the same arguments will apply to the integral (|3. 10b[) ] . 
The phase <p(t) [see Eq. (|2.6j) ] can be written as 



0(t) =0J [t + (f) m (t)]+<f> s (t), 



where 



0m (t) : = 



n ■ r d (i) 



(6.1) 



3.2a) 



k=l 



(k + iy. 



n -r d (i)<A t k 

— - — 2^ u *w (6 - 2b) 



k=l 



The functions 4> m (t) and <fi s {t) do not depend on fre- 
quency ujq. We can write the integral (|3.10a[) as 



F„ 



x(t) a(t) e ^ s(t) exp { - kj [t + <p m (t)] } dt. 



(6.3) 

Next we introduce a new time variable so called 
barycentric time [l], |33| , 



t h (t) :=t + c/> m (t). 



(6.4) 



In this new time coordinate the integral (|6.3I) is approx- 
imately given by (see Ref. [l|, Sec. Ill D) 



F a £* 



x[t(t h )}a[t(t b )}( 



icj) e [t(tb)]g— iwotb 



dt h 



(6.5) 



This integral is a Fourier transform of the data £c[t(tb)] 
multiplied by the function a[i(tb)] exp[— i0 s [t(tb)]]. For 
discrete data x{t) the integral (|6.5|) can be converted to 
a discrete Fourier transform which can be evaluated by 
the FFT algorithm. 

Thus to convert the integral (|6.3p into a Fourier trans- 
form we need to resample the function x(t) a(t) e -1 ^*) 
according to Eq. (|6.4p . We consider two numerical inter- 
polation methods in order to obtain the resampled func- 
tion. The first method is the nearest neighbor interpola- 
tion also called the stroboscopic resampling. We assume 
that the original data is a time series Xk (k = 1, . . . , N), 
sampled at uniform intervals. In this method we obtain 
the value of the time series Xk at barycentric time t b by 
taking the value yk such that fco is the nearest integer 
to %. We have illustrated the method in Fig. [3] 

The second method has two steps. The first step con- 
sists of obtaining a more finely sampled time series and 
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FIG. 4: Comparison of the two interpolation methods and the 
perfect matched-filtering. We see that the two-step interpo- 
lation method that uses Fourier and spline interpolation very 
accurately reproduces the perfect matched filter. 



-1 01234567 
Time 



FIG. 3: Illustration of the nearest neighbor resampling 
method. The top panel shows the uniformly sampled orig- 
inal data and the points of the barycentric time. The middle 
panel shows the interpolations of the original time series at 
the points of the barycentric time, obtained by the nearest 
neighbor method. The bottom panel shows the resampled 
time series which is the uniformly sampled barycentric time 
series form the middle panel. 



the second step consists of interpolating the upsampled 
time series to the barycentric time using splines (34|. To 
perform the first step we use an interpolation method 
based on the Fourier transform. We take the Fourier 
transform of the original time series, pad the Fourier 
transform series with an appropriate amount of zeros 
and then transform it back to the time domain by in- 
verse Fourier transform. The Fourier transforms are per- 
formed using the FFT algorithm. We thus obtain an 
interpolated time series with points inserted between the 
original points. If we have a time series with N points 
and pad its discrete Fourier transform with N zeros, by 
inverse transform we obtain a 2A r -point time series. The 
second step consists in applying splines to interpolate the 
upsampled time series to the barycentric time for num- 
ber of points equal to the number of original data points. 
Thus if the original time series contains N points the final 



interpolated time series contains also N points. 

We have compared the performance of the two interpo- 
lation methods and we have also compared these meth- 
ods with an exact matched filter. To carry out the com- 
parison we have used noise-free waveforms given in Eqs. 
(|2.1[) - (l2.6[) with one spin-down parameter. We have cal- 
culated the ^-statistic using the two interpolation meth- 
ods and exact matched-filtering method. In the matched 
filtering method we have assumed that we know the fre- 
quency of the signal and thus the Doppler modulation 
due to the motion of the detector. However, we have 
used FFT to calculate the F-statistic for the whole fre- 
quency band. The results are shown in Fig. 2J 

We have performed a Monte Carlo simulation consist- 
ing of 1000 trials and we have found that the rms er- 
ror divided by maximum of the F-statistic in the second 
method was 0.1% whereas in the first, fastest method 
it was 5%. The nearest neighbor interpolation leads to 
a greater signal-to-noise ratio loss than spline interpola- 
tion and also, very importantly, to elevated sidelobes of 
the F-statistic. In the presence of noise this can lead to 
a loss of the parameter estimation accuracy if the noise 
elevates the sidelobes above the main maximum. The 
stroboscopic resampling is much faster than the second 
two step method however the second method is much 
more accurate than the first. 



B. FFT interpolation 

Using the FFT algorithm we can efficiently calcu- 
late the discrete Fourier transform (DFT) X{k) {k = 
1, . . . , N) of a time series xi (£ = 1, . . . , N). We recall 
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that X(k) is given by the following expression 

N 

X(k) = J2w 2ni(i ~ m ~ 1)/N i k = l,...,N. (6.6) 

1=1 

The frequencies (fc — 1)/N are called Fourier frequencies 
and DFT components calculated at Fourier frequencies 
are called Fourier bins. When the true frequency of a 
monochromatic signal does not coincide with one of the 
Fourier frequencies, the use of the FFT algorithm to eval- 
uate the sum (|6.6p leads a certain loss of signal-to-noise 
ratio. The greatest loss equal around 36.3% is when the 
true frequency is half way between the Fourier frequen- 
cies. 

One way to improve this situation is to pad the time 
series of N data points with N zeros. This leads to DFT 
evaluated at twice as many points as the DFT of the 
original time series and the signal-to-noise loss is only 
9.97%. However this procedure leads to evaluating twice 
as long FFT as the original ones and thus increases the 
computational time by more than a factor of two. 

There exists an approximate interpolation procedure 
proposed by pulsar astronomers (see Chapter 7.3.3 in 
Ref. [2l|), in which the DFT component in the middle 
of two Fourier frequencies is approximated by 

X(k + l/2) = [X(k + l)- X(k)}/V2. (6.7) 

This interpolation method is called interbinning. One 
can show (see Fig. 7.3 in Ref. [HI]) that the interpolation 
based on Eq. (|6.7|) leads to maximum loss of signal-to- 
noise ratio of 13%. 



C. Finding the maximum of the ^-statistic 
accurately 

As we calculate the ^-statistic on a discrete grid in the 
parameter space and as the maximum of the ^-statistic 
does not in general coincide with a node of the grid, we 
always lose some signal-to-noise ratio and parameter es- 
timation accuracy. To improve this we use a nonlinear 
optimization routine to find an improved maximum of 
our statistic. 

The search for the maximum of the ^-statistic is per- 
formed in two steps. First we find the maximum of T 
over the discrete grid in the parameter space and then 
the parameters obtained in this coarse search we input 
as initial values to some hill-climbing optimization rou- 
tine to find an improved maximum. The second step is 
called a fine search. As our maximum finding routine we 
can use a direct search that does not require calculations 
of derivatives of J 7 , namely the Nelder-Mead algorithm 
or simplex search algorithm 36]. The Nelder-Mead al- 
gorithm is illustrated in Fig. [5] for two-dimensional case 
where simplices are triangles. 



Nelder-Mead algorithm 
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FIG. 5: An illustration of the Nelder-Mead algorithm in two 
dimensions. We first evaluate the function to be minimized on 
an initial triangle. At each step of the search, a new point in 
or near the current triangle is generated. The function value 
at the new point is compared with the function's values at 
the vertices of the triangle and, usually, one of the vertices is 
replaced by the new point, giving a new triangle. This step 
is repeated until the diameter of the simplex is less than the 
specified tolerance. In the method the triangle adapts itself 
to the local landscape, elongating down long inclined planes, 
changing direction on encountering a valley at an angle, and 
contracting in the neighborhood of a minimum. 



VII. MONTE CARLO SIMULATIONS 

We have implemented the data analysis methods pre- 
sented in the previous sections in a computer code and 
we have performed a number of Monte Carlo simulations 
to test how accurately we can estimate signal's parame- 
ters for different values of the signal-to-noise ratio (SNR) 
and how the rms errors of our ML estimators compare to 
the Cramer-Rao lower bound. For unbiased estimators 
the Cramer-Rao lower bound on variances of the esti- 
mators is given by the diagonal elements of the inverse 
of the Fisher matrix. The ML estimators are asymp- 
totically (i.e. when SNR tends to infinity) unbiased and 
with variances approaching the diagonal elements of the 
inverse of the Fisher matrix. Our Monte Carlo simula- 
tions consisted of generating the data that were the sum 
of the white noise and the gravitational-wave signal pre- 
sented in Sec. |TT] and using the filtering procedure from 
Sec. Mil to detect the signal and estimate its parameters. 
We have taken the observation time equal to exactly 2 
sidereal days. For this observation time it is enough to 
take only the first spin-down parameter in the templates 
(seeRefs. @,[|). To obtain the Cramer- Rao lower bound 
we have calculated the 8x8 Fisher matrix [defined in Eq. 
(|3.17[) ] for the signal model given in Sec.|TT]with one spin 
down parameter included. 
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We have made our simulation to mimic the analysis 
of the VIRGO interferometer VSR1 data. We have thus 
assumed that data comes from an interferometric detec- 
tor located at the position of the VIRGO detector and 
we have used the detector's ephemerides corresponding 
to data taken sometime in the year 2007. We have gen- 
erated data in a narrow frequency band of 50 mHz with 
the lower edge of the band equal to 435.1875 Hz. This re- 
sulted in a short time series of 17233 data points greatly 
reducing the CPU time needed to perform the simula- 
tions. For all the simulations we have chosen the same 
parameters of the gravitational- wave signal except for the 
constant amplitude ho that we scaled to obtain data with 
a chosen SNR [defined in Eq. ()3.16|) ] . We have found the 
grid point p g nearest to the true position of the signal 
in the parameter space and we have calculated the F- 
statistic on a small grid around the point p g . The size 
of the grid was ±2 grid points from the point p g in the 
direction of the parameters Lb, a±, ct2 [see Eqs. Q3.12)) ]. 
For each set of these three parameters we have evaluated 
the /"-statistic for the whole 50 mHz band. 

In all the simulations we have used the constrained grid 
constructed in Sec. IIVI and we have chosen the detection 
threshold for the F-statistic equal to 10. This low thresh- 
old ensured that the probability of detection was nearly 
one and none of the signals was missed. In each simula- 
tion the estimation of signal's parameters was performed 
in two steps. The first step, called the coarse search, was 
the calculation of the .F-statistic on the grid described in 
Sec. IIVI We have registered all the threshold crossings 
and we have taken the coarse estimates of the signal's 
parameters as the parameters of the grid point for which 
the F-statistic was maximal. In the second step, called 
the fine search, we have found the maximum of the F- 
statistic for each signal registered in the first step using 
the Nelder-Mead algorithm with the initial values equal 
to the coarse estimates of the parameters. In this step 
the F-statistic was calculated [by means of Eqs. (|3.9[) - 
(I3.10[) ] exactly without any approximations. For each 
value of the SNR we have repeated the simulation 1000 
times with different realizations of the white Gaussian 
noise. 

In the first simulation we have employed a fine grid 
with the minimal match MM = \/0.9. We have used the 
accurate two-step resampling procedure described in Sec. 
IVI Al and before applying the FFT we have padded the 
data with 15535 zeros resulting in a time series of 2 15 
points. This ensured an exact interpolation of the DFT 
between the Fourier frequencies and the fastest imple- 
mentation of the FFT algorithm (because the number of 
data points was a power of 2). The results of this sim- 
ulation are presented in Fig. [6l where we have depicted 
standard deviations and biases of the intrinsic parame- 
ters of the signal: frequency, spin down, declination, and 
right ascension, as functions of the SNR. The results of 
simulation of 1000 runs are marked by the circles. Addi- 



tionally in the left panels we have shown the Cramer- Rao 
bounds for the standard deviations calculated from the 
inverse of the Fisher information matrix. We see that for 
the SNRs greater than ~9 the standard deviations ob- 
tained from simulation are very close to the Cramer-Rao 
bounds and the simulated biases are small fractions of a 
percent of the true values. 

In the second simulation we have tested how accuracy 
of the parameter estimation is affected by various options 
of the algorithms described in the previous section. All 
these options aim at speeding up computations. In the 
simulation we were changing the thickness of the grid, we 
were comparing the two-step spline resampling with the 
nearest neighbor resampling, and we were also testing 
the interbinning interpolation [defined in Eq. (16. 7[) ]. In 
all the runs we have used a thicker grid corresponding to 
the minimal match MM — v / 3/2. We have studied three 
specific cases: (i) zero padding and spline interpolation; 
(ii) interbinning and spline interpolation; (iii) interbin- 
ning and the nearest neighbor interpolation. The results 
of the simulation are presented in Fig. [7] From compar- 
ison of the three cases it follows that with a coarser grid 
and with the use of different approximations the rms er- 
rors of the parameter estimators are greater but still at 
a reasonable level. In the range of the SNRs from —9 to 
~15, the standard deviations obtained from simulation 
are twice as large as the Cramer-Rao bounds, when we 
use spline resampling. For the nearest neighbor interpo- 
lation, for some parameters, the simulated rms error is 
twice as large as the Cramer- Rao bound for the SNRs up 
to -30. 

In the third simulation we have studied another three 
specific cases: (i) interbinning and spline interpolation, 
with a coarse grid of MM = v/3/2; (ii) interb inning with 
the nearest neighbor interpolation, with a fine grid of 
MM = ^0.9; (iii) interbinning with the nearest neighbor 
interpolation, with a coarse grid of MM = -\/3/2. The re- 
sults of the simulation are presented in Fig. |SJ This sim- 
ulation shows that with a sufficiently fine grid even the 
use of the least accurate (but the fastest) resampling by 
the nearest neighbor interpolation leads to the rms errors 
of the intrinsic parameters very close to the Cramer-Rao 
bounds. 



Acknowledgments 

The contributions of K. M. Borkowski, P. Jaranowski, 
A. Krolak, and M. Pictka were supported in part by the 
MNiSzW grants nos. 1 P03B 029 27 and N N203 387237. 
A.K. would like to acknowledge hospitality of the Max 
Planck Institute for Gravitational Physics in Hannover, 
Germany, where part of this work was done. We would 
also like to thank Holger Pletsch for discussions and help- 
ful remarks. 



17 




FIG. 6: Standard deviations (left panels) and biases (right panels) of the ML estimators of the intrinsic parameters as functions 
of the SNR computed with the grid of MM = ^0.9. The plots are made, from the top, for the following parameters: frequency, 
spin down, declination, and right ascension. Each circle corresponds to a standard deviation or a bias computed from simulation 
of 1000 runs. The continuous lines in the left panels are the Cramer-Rao bounds for standard deviations. 



18 




FIG. 7: Standard deviations (left panels) and biases (right panels) of the ML estimators of the intrinsic parameters as functions 
of the SNR computed with the grid of MM = for the following three cases: (i) zero padding and spline interpolation 

(diamonds); (ii) interbinning and spline interpolation (stars); (iii) interbinning and the nearest neighbor interpolation (circles). 
The continuous lines in the left panels are the Cramer-Rao bounds. 
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FIG. 8: Standard deviations (left panels) and biases (right panels) of the ML estimators of the intrinsic parameters as functions 
of the SNR computed for the following three cases: (i) interbinning and spline interpolation, with the coarse grid of MM = s/3/2 
(diamonds); (ii) interbinning with the nearest neighbor interpolation, with the fine grid of MM = (stars); (iii) interbinning 

with the nearest neighbor interpolation, with the coarse grid of MM = \/3/2 (circles). The continuous lines in the left panels 
are the Cramer-Rao bounds. 
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